Population Parameter Estimation

Estimation of pediatric atorvastatin PK given informative priors based on adults




Specification of prior distributions for population means and inter-individual covariance

Let’s revise the priors to be more consistent with current practice. Specifically,

  • Let’s use weakly informative priors to guard against implausible values and help improve sampling performance.
    • Update variance of CL/F, V2/F, and Q/F to 2.
    • Add variance of \(F_{chew}\) and SD to 1.
  • Let’s use a LKJ prior for IIV correlation matrix and lognormal priors for diag(OMEGA).
    • OLKJDF = 1 and OVARF=1.

Setup R environment

rm(list = ls())
gc()

modelName <- "atorv1"
scriptName <- paste(modelName, "Rmd", sep = ".")
fitModel <- FALSE

## Relative paths assuming the working directory is the script directory
## containing this script
scriptDir <- getwd()
projectDir <- dirname(scriptDir)
dataDir <- file.path(projectDir, "data")
modelDir <- file.path(projectDir, "model")
outDir <- file.path(modelDir, modelName)
figDir <- file.path(projectDir, "deliv", "figure", modelName)
tabDir <- file.path(projectDir, "deliv", "table", modelName)
invisible(dir.create(figDir, recursive = TRUE))
## Warning in dir.create(figDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/figure/atorv1' already exists
invisible(dir.create(tabDir, recursive = TRUE))
## Warning in dir.create(tabDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/table/atorv1' already exists
.libPaths("lib")

library(metrumrg)
## Loading required package: reshape
## Loading required package: lattice
## Loading required package: XML
## Loading required package: MASS
## Loading required package: grid
## metrumrg 5.57
## enter "?metrumrg" for help
suppressMessages(library(rstan))
suppressMessages(library(bayesplot))
suppressMessages(library(tidyverse))
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(parallel)
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(qapply)
## Loading required package: rlecuyer
## 
## Attaching package: 'qapply'
## The following object is masked from 'package:metrumrg':
## 
##     qstat
library(PKPDmisc)

knitr::opts_chunk$set(
  comment = '.', 
  fig.height = 5, 
  fig.width = 9,
  eval.after = 'fig.cap',
  message = FALSE,
  fig.path = file.path(figDir, paste(modelName, "_", sep = ""))
)

## Go back to default ggplot2 theme that was overridden by bayesplot
theme_set(theme_gray())

source(file.path(scriptDir, "tools", "stanTools2.R"))
source(file.path(scriptDir, "tools", "functions.R"))

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

set.seed(11191951) 
set.seed(10271998) 

The NONMEM model stub

modelFile <- file.path(modelDir, paste(modelName, "stub.ctl", sep = ""))
modelText <- readLines(modelFile)
cat(modelText, sep = "\n")
. $PROBLEM RUN# atorv1
. $INPUT C ID TIME AMT DV CMT DOSE EVID MDV ADDL II SS FORM WT PER BLQ
. ;$INPUT C ID SID=DROP STDY=DROP TIME DV DVR=DROP DVC=DROP AMT AMTO=DROP CMT TYPE=DROP
. ;       DOSE DOSN=DROP NTIM=DROP EVID MDV ADDL II SS FORM WT SEX=DROP AGE
. ;       HT=DROP RACE=DROP PER BLQ TRCD=DROP ANLT=DROP TXT=DROP 
. ;$DATA ../../../data/atorv1172atvb.csv IGNORE=(C='C', BLQ.EQ.1)
. $DATA ../../../data/atorvWrkShop.csv IGNORE=(C='C', BLQ.EQ.1)
. 
. $ABBR DECLARE INTEGER FIRST_WRITE_PAR
. $ABBR DECLARE INTEGER FIRST_WRITE_IPAR
. 
. $SUBROUTINES ADVAN4 TRANS4 OTHER = ../../extrasend.f90
. 
. $PK
. ; Request extra information for Bayesian analysis.
. ; An extra call will then be made for accepted samples 
. include '/opt/NONMEM/nm74/util/nonmem_reserved_general'
. BAYES_EXTRA_REQUEST=1
. 
. nThin = THETA(8)
. 
. VWT = LOG(WT/70) ; normalized to 70 kg adult
. MU_1 = THETA(1) + VWT * 0.75 ; CL
. MU_2 = THETA(2) + VWT        ; V2
. MU_3 = THETA(3) + VWT * 0.75 ; Q
. MU_4 = THETA(4) + VWT        ; V3
. MU_5 = THETA(5)              ; ka
. MU_6 = THETA(6)              ; F1chew
. MU_7 = THETA(7)              ; SD
. 
. " CALL EXTRASEND()
. 
. ;ATORVASTATIN
. 
. CL = EXP(MU_1 + ETA(1)) ; atv cl 
. V2 = EXP(MU_2 + ETA(2))
. Q = EXP(MU_3 + ETA(3))
. V3 = EXP(MU_4 + ETA(4))
. ;deltaKA = EXP(MU_5 + ETA(5))
. KA = EXP(MU_5 + ETA(5))
. F1chew = EXP(MU_6 + ETA(6))
. SD = EXP(MU_7 + ETA(7))
. 
. ; Constrain KA > lambda1 (typical value) to enhance identifiability
. ;T1 = CL/V2
. ;T23 = Q/V2
. ;T32 = Q/V3
. ;L1 = ((T1+T23+T32)+SQRT((T1+T23+T32)**2-4*T1*T32))/2
. ;L2 = ((T1+T23+T32)-SQRT((T1+T23+T32)**2-4*T1*T32))/2
. 
. ;KA = deltaKA + L2
. S2 = V2/1000              ; DOSE IN uM & CONC IN nM/L
. F1 = 1.0             ; TABLET AS REFERENCE
. IF(FORM.EQ.2) F1 = F1chew  ; CHEWABLE F1
. 
. IF(BAYES_EXTRA==1 .AND. NIREC==1 .AND. NDREC==1 .AND. &
.    ITER_REPORT>0 .AND. &
.    ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN 
.   IF(FIRST_WRITE_PAR==0) THEN
.     " OPEN(unit=52,FILE='./par.txt') 
.     FIRST_WRITE_PAR=1
.   ENDIF
.   " WRITE(52,'(I12,1X,20(1X,1PG19.10E3))') ITER_REPORT, &  
.   " OMEGA(1,1), OMEGA(2,1), OMEGA(2,2), &
.   " THETA(1), THETA(2), THETA(3), THETA(4), THETA(5), THETA(6), THETA(7)
. ENDIF
. IF(BAYES_EXTRA==1 .AND. NDREC==1 .AND. ITER_REPORT>0 .AND. &
.   ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN 
.   IF(FIRST_WRITE_IPAR==0) THEN
.     " OPEN(unit=50,FILE='./ipar'//TFI(PNM_NODE_NUMBER)//'.txt') 
.     FIRST_WRITE_IPAR=1
.   ENDIF
.   " WRITE(50,'(I12,1X,F14.0,10(1X,1PG19.10E3))') ITER_REPORT, ID, &
.   " ETA(1), ETA(2)
. ENDIF
. 
. 
. $ERROR
. 
. LOQ = 0.45          ; nM/L
. IPRED = F
. DUM = (LOQ-IPRED)/(SD*IPRED)
. CUMD=PHI(DUM)
. IF(BLQ.EQ.0)THEN
.   F_FLAG=0
.   Y = F*(1+SD*ERR(1)) ; ATV
. ELSE
.   F_FLAG=1
.   Y=CUMD
. ENDIF
. 
. ; Code below into chains

Assemble initial estimates, priors and table specs

## create initial estimates
geninit <- function(){
    paste(c(
    "; Initial values of THETA",
    "$THETA",
    paste(rnorm(1, log(1000), 0.5), "     ; TOTAL CL ATV = THETA(1)"),
    paste(rnorm(1, log(600), 0.5), "     ; V2 = THETA(2)"),
    paste(rnorm(1, log(130), 0.5), "   ; Q = THETA(3)"),
    paste(rnorm(1, log(1500), 0.5), "   ; V3 = THETA(4) "),
    paste(rnorm(1, log(1.2), 0.5), "   ; KA= THETA(5) "),
    paste(rnorm(1, log(1), 0.25), "   ; THETA(6) relative F1 CHEWABLE = THETA(5) "),
    paste(log(0.2), " ; THETA(7) - PRO RES ERROR"),
    "$OMEGA BLOCK(2) ;INITIAL values of OMEGAs",
    "0.1003           ; cl",
    ".06034  0.8477     ; v2",
    "$OMEGA 0 FIX",
    "$OMEGA 0 FIX",
    "$OMEGA 0 FIX",
    "$OMEGA 0 FIX",
    "$OMEGA 0 FIX",
    "$SIGMA  ;Initial value of SIGMA",
    "1 FIX; pro"),
    sep = "\n")
}

## Set parameters of the prior distribution
priors <- paste(c("$PRIOR NWPRI",
                  "$THETAP          ; Prior information of THETAS",
                  paste("(", 5.65, " FIX)      ;  THETA(1) TOTAL CL ATV"),
                  paste("(", 5.60, " FIX)      ;  THETA(2) V2"),
                  paste("(", 4.29, " FIX)      ;  THETA(3) Q"),
                  paste("(", 7.07, " FIX)      ;  THETA(4) V3"),
                  paste("(", -1.33, " FIX)      ; THETA(5) KA"),
                  paste("(", log(1), " FIX)      ; THETA(6) F1 chewable"),
                  paste("(", log(0.2), " FIX)      ; THETA(7) log(sigma)"),
                  "$THETAPV BLOCK(7)     ;  variances for priors on THETAS (var-cov)",
                  "2 FIX ; CL weakly informative",
                  "0.00 2  ; V2 weakly informative",
                  "0.00 0.00 2  ; Q weakly informative",
                  "0.00 0.00 0.00 .0156   ; V3 SE**2 informative",
                  "0.00 0.00 0.00 0.00 .00684    ; KA SE**2 informative",
                  "0.00 0.00 0.00 0.00 0.00 1    ; F1 chewable weakly informative",
                  "0.00 0.00 0.00 0.00 0.00 0.00 1    ; log(sigma) weakly informative",
## Parameters LKJ / lognormal prior; off-diagonal and OMEGAP ignored
                  "$OMEGAP BLOCK(2)   ; prior information for OMEGA",
                  ".1003   FIXED         ; cl",
                  ".06034    .8477       ; v2" ,
                  "$OMEGAPD (4.0, FIXED)     ; df for OMEGA prior"
## Parameters for inverse Wishart prior
##                  "$OMEGAP BLOCK(2)   ; prior information for OMEGA",
##                  ".1003   FIXED         ; cl",
##                  ".06034    .8477       ; v2" ,
##                  "$OMEGAPD (4.0, FIXED)     ; df for OMEGA prior"
),
sep = "\n")

## Table statements
tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE WT",
                  "NOPRINT ONEHEADER FILE=./.tab",
                  "$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
                  "NOPRINT ONEHEADER FILE=./par.tab"),
                sep = "\n")
##tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE AGE WT",
##                  "NOPRINT ONEHEADER FILE=./.tab",
##                  "$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
##                  "NOPRINT ONEHEADER FILE=./par.tab"),
##                sep = "\n")

Run NONMEM

nChains <- 4
nPost <- 250 ## Number of post-burn-in samples per chain after thinning
nBurn <- 250 ## Number of burn-in samples per chain after thinning
nThin <- 1

seed = sample(10000, size = nChains)
if(fitModel){
  if(!file.exists(file.path(modelDir, modelName))) 
    dir.create(file.path(modelDir, modelName))
  runChains <- mclapply(1:nChains, runChain,
                      modelName = modelName, modelDir = modelDir, 
                      priors = priors, tables = tables,
                      nPost = nPost, nBurn = nBurn, nThin = nThin, 
                      seed = seed, print = 100,
                      ## LKJ prior for IIV correlation matrix with oarameter = OLKJDF
                      OLKJDF = 1, 
                      ## lognormal prior for IIV SDs
                      ## log(sqrt(Omega(i))) ~ Normal(log(sqrt(OmegaPrior(i))), 1/OVARF)
                      OVARF = 1, AUTO = 2, 
                      NUTS_DELTA = 0.8,
                      grid = FALSE,
                      method = "NUTS", 
##                      pe = "orte 16",
                      mode = "nonpara")
  # NONR(run = paste(modelName, ".", 1, sep = ""),
  #        command = "/opt/NONMEM/nm74/nmqual/autolog.pl",
  #        project = file.path(modelDir, modelName),
  #        grid = FALSE,
  #        wait = FALSE,
  #        diag = FALSE,
  #        fdata = FALSE,
  #        purge = FALSE,
  #        checkrunno = TRUE)
  
  }

read in data file

dataFile <- "atorvWrkShop.csv"
nmData <- read.csv(file.path(dataDir, dataFile), as.is = TRUE)

Get population parameters

chains <- 1:nChains
##chains <- c(1, 4)

nChains2 <- length(chains)

## Read in posterior distributions by chain and add a column called chain
popParameters <- map(chains, function(thisChain){
  param <- data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""), 
                                       "par.txt"), data.table = FALSE) 
  names(param) <- c("iteration", "OMEGA11", "OMEGA21", "OMEGA22",
                    "THETA1", "THETA2", "THETA3", "THETA4", "THETA5", 
                    "THETA6", "THETA7")
  param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
  mutate(sample = 1:n())

## Calculate more interpretable paraneters
popParameters2 <- popParameters %>%
  mutate(CLHat = exp(THETA1),
         V2Hat = exp(THETA2),
         QHat = exp(THETA3),
         V3Hat = exp(THETA4),
##         k10 = CLHat / V2Hat,
##         k12 = QHat / V2Hat,
##         k21 = QHat / V3Hat,
##         lambda = ((k10 + k12 + k21) - 
##           sqrt((k10 + k12 + k21)^2 - 4 * k10 * k21)) / 2,
##         KAHat = exp(THETA4) + lambda,
##         KAHat = exp(THETA4 + lambda),
         KAHat = exp(THETA5),
         F1Chewable = exp(THETA6),
         sigma1 = exp(THETA7)) %>%
  select(starts_with("OMEGA"),
         CLHat, V2Hat, QHat, V3Hat, KAHat, F1Chewable, sigma1, chain, iteration, sample)

## Convert to 3-D array with dims = {iterations, chains, parameters}. 
popParArray <- array(double(nPost * nChains2 * (ncol(popParameters2) - 3)), 
                     dim = c(nPost, nChains2, ncol(popParameters2) - 3), 
                     dimnames =  list(NULL, NULL, setdiff(names(popParameters2), c("chain", "iteration", "sample"))))
for(iChain in 1:nChains2){
  popParArray[,iChain,] <- popParameters2 %>%
    filter(chain == chains[iChain],
           iteration > 0) %>%
    select(-chain, -iteration, -sample) %>%
    as.matrix
}

Get individual parameters

indParameters = map(chains, function(thisChain){
  files <- list.files(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = "")))
  iparFiles <- files[grepl("^ipar", files)]
  param <- map(iparFiles, function(thisFile){
    data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""), 
                                thisFile), data.table = FALSE)
  }) %>% 
    bind_rows 
  names(param) <- c("iteration", "ID", "ETA1", "ETA2")
  param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
  arrange(chain, iteration, ID) %>%
  left_join(popParameters)

MCMC diagnostics and posterior distributions of parameters

options(bayesplot.base_size = 12,
        bayesplot.base_family = "sans")
color_scheme_set(scheme = "brightblue")
myTheme <- theme(text = element_text(size = 12), 
                 axis.text = element_text(size = 12))

plot_mcmcHistory <- mcmcHistory(popParArray,
                                 nParPerPage = 5, myTheme = myTheme)
plot_mcmcDensityByChain <- mcmcDensity(popParArray, 
                                        nParPerPage = 16, byChain = TRUE, 
                                        myTheme = theme(text = element_text(size = 12), 
                                                        axis.text = element_text(size = 10)))
plot_mcmcDensity <- mcmcDensity(popParArray, nParPerPage = 16, 
                                 myTheme = theme(text = element_text(size = 12), 
                                                 axis.text = element_text(size = 10)))

plot_mcmcHistory
. [[1]]
Pediatric atorvastatin PK example: MCMC trace plots.

Pediatric atorvastatin PK example: MCMC trace plots.

. 
. [[2]]
Pediatric atorvastatin PK example: MCMC trace plots.

Pediatric atorvastatin PK example: MCMC trace plots.

plot_mcmcDensityByChain
. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities by chain.

Pediatric atorvastatin PK example: Univariate marginal densities by chain.

plot_mcmcDensity
. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities.

Pediatric atorvastatin PK example: Univariate marginal densities.

mcmc_pairs(popParArray[,,!grepl("^OMEGA", dimnames(popParArray)[[3]])])
Pediatric atorvastatin PK example: Bivariate marginal densities.

Pediatric atorvastatin PK example: Bivariate marginal densities.

mcmc_pairs(popParArray[,,grepl("^OMEGA", dimnames(popParArray)[[3]])])
Pediatric atorvastatin PK example: Bivariate marginal densities.

Pediatric atorvastatin PK example: Bivariate marginal densities.

caption <- paste("Pediatric atorvastatin PK example:",
      c(rep("MCMC trace plots.", length(plot_mcmcHistory)),
        rep("Univariate marginal densities by chain.",
            length(plot_mcmcDensityByChain)),
        rep("Univariate marginal densities.",
            length(plot_mcmcDensityByChain)),
        rep("Bivariate marginal densities.", 2)))

ptable <- monitor(popParArray, 
                  warmup = 0, print = FALSE) %>% 
  as.matrix %>%
  formatC(3) %>%
  as.data.frame
write.csv(ptable, file = file.path(tabDir, paste(modelName, 
                                                 "ParameterTable.csv", sep = "")))

ptable %>% 
  rename(SEmean = se_mean, SD = sd, pct2.5 = "2.5%", pct25 = "25%", median = "50%",
         pct75 = "75%", pct97.5 = "97.5%", Neff = "n_eff") %>%
  mutate(parameter = rownames(.), "95% CI" = paste("(", pct2.5, ", ", pct97.5, ")", 
                                                   sep = "")) %>%
  select(parameter, mean, SD, median, "95% CI", Neff, Rhat) %>%
  kable(caption = "Summary of model parameter estimates.") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Summary of model parameter estimates.
parameter mean SD median 95% CI Neff Rhat
OMEGA11 0.416 0.171 0.419 (0.195, 0.769) 3 1.56
OMEGA21 0.0898 0.128 0.125 (-0.213, 0.305) 37 1.11
OMEGA22 294 508 1.04 (0.48, 1.17e+03) 2 1.53
CLHat 590 151 618 ( 373, 842) 3 1.57
V2Hat 2.44e+03 2.36e+03 1.17e+03 ( 667, 6.52e+03) 2 1.53
QHat 394 213 293 ( 188, 755) 2 1.52
V3Hat 1.6e+03 200 1.6e+03 (1.37e+03, 2.01e+03) 4 1.33
KAHat 0.27 0.0396 0.281 (0.209, 0.329) 2 1.53
F1Chewable 0.896 0.179 0.87 (0.675, 1.21) 374 1.55
sigma1 0.258 0.0494 0.235 (0.208, 0.341) 2 1.53

Construct design data frame

## For now just use the NONMEM data set

design <- nmData %>%
  filter(C == ".")

mrgsolve model code for posterior predictions

library(mrgsolve)

indCode =
  '
$PARAM 
THETA1 = 5.65
THETA2 = 5.60
THETA3 = 4.29
THETA4 = -1.33
THETA5 = 7.07
THETA6 = 1
THETA7 = -3.2189
WT = 70
FORM = 1
ETA1 = 0
ETA2 = 0

$MAIN

double VWT = log(WT/70); 
double CL = exp(THETA1 + (VWT * 0.75) + ETA1 + ETA(1)); 
double V2 = exp(THETA2 + VWT + ETA2 + ETA(2));
double Q = exp(THETA3 + (VWT * 0.75));
double V3 = exp(THETA5 + VWT);

//double T1 = CL/V2;
//double T23 = Q/V2;
//double T32 = Q/V3;
//double L2 = ((T1+T23+T32)-sqrt(pow(T1+T23+T32, 2)-4*T1*T32))/2;

double KA = exp(THETA4);
F_GUT = 1.0;             // TABLET AS REFERENCE
if(FORM == 2) F_GUT = exp(THETA6);  // CHEWABLE F1

$PKMODEL cmt = "GUT CENT PERIPH", depot = TRUE

$SIGMA 1

$OMEGA @block
0 0 0

$TABLE
double cPredHat = 1000 * (CENT / V2);
double cPred = cPredHat * (1 + exp(THETA7) * EPS(1));

$CAPTURE cPred
'

## Compile model
atvMod <- mcode("ATVmodel", indCode)

if(FALSE){
  ## Test runs
  out <- atvMod %>%
    data_set(design) %>%
    carry_out(DOSE) %>%
    mrgsim
  
  plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
  
  out <- atvMod %>%
    data_set(design) %>%
    param(popParameters2[1,] %>%
            rename(CL = CLHat,
                   V2 = V2Hat,
                   Q = QHat,
                   V3 = V3Hat,
                   KA = KAHat,
                   F1 = F1Chewable)) %>%
    carry_out(DOSE) %>%
    mrgsim
  
  plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
}

Individual predictions

samples <- unique(indParameters$sample)
samples <- 1:8
indSim1 <- mclapply(samples,
                    function(thisSample){
                      atvMod %>%
                        data_set(design) %>%
                        idata_set(indParameters %>% 
                                filter(sample == thisSample)) %>%
                        carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
                        mrgsim %>%
                        as.tbl
                    }) %>%
  bind_rows

Population predictions

samples <- unique(popParameters$sample)
popSim1 <- mclapply(samples,
                    function(thisSample){
                      atvMod %>%
                        data_set(design) %>%
                        param(popParameters %>% 
                                filter(sample == thisSample)) %>%
                        omat(matrix(with(popParameters %>% 
                                filter(sample == thisSample), 
                                c(OMEGA11, OMEGA21, OMEGA21, OMEGA22)), ncol = 2)) %>%
                        carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
                        mrgsim %>%
                        as.tbl
                    }) %>%
  bind_rows

Posterior predictive distributions

## Combined plots

indPred <- indSim1 %>%
  group_by(ID, TIME, EVID, DOSE, FORM) %>%
  summarize(lbInd = quantile(cPred, probs = 0.05, na.rm = TRUE),
            medianInd = quantile(cPred, probs = 0.5, na.rm = TRUE),
            ubInd = quantile(cPred, probs = 0.95, na.rm = TRUE))

popPred <- popSim1 %>%
  group_by(ID, TIME, EVID, DOSE, FORM) %>%
  summarize(lbPop = quantile(cPred, probs = 0.05, na.rm = TRUE),
            medianPop = quantile(cPred, probs = 0.5, na.rm = TRUE),
            ubPop = quantile(cPred, probs = 0.95, na.rm = TRUE))

predAll <- indPred %>%
  left_join(popPred) %>%
  left_join(design) %>%
  mutate(DV = suppressWarnings(as.numeric(DV)))

forms <- sort(unique(predAll$FORM))
plot_PPCPK <- lapply(forms,
                     function(thisForm){
  p1 <- ggplot(predAll %>% filter(FORM == thisForm), 
               aes(x = TIME, y = DV))
  p1 <- p1 + 
    geom_line(aes(x = TIME, y = medianPop, 
                 color = "population")) +
    geom_ribbon(aes(ymin = lbPop, ymax = ubPop, 
                   fill = "population"), 
               alpha = 0.25) +
    geom_line(aes(x = TIME, y = medianInd, 
                  color = "individual")) +
    geom_ribbon(aes(ymin = lbInd, ymax = ubInd, 
                    fill = "individual"), 
                alpha = 0.25) +
    scale_color_brewer(name  ="prediction",
                       breaks=c("individual", "population"),
                       palette = "Set1") +
    scale_fill_brewer(name  ="prediction",
                      breaks=c("individual", "population"),
                      palette = "Set1")
  p1 + geom_point() +
    labs(title = paste("formulation =", c("tablet", "chewable")[thisForm]),
         x = "time (h)",
         y = "atorvastatin plasma concentration (nM)") +
    theme(text = element_text(size = 12),
          axis.text = element_text(size = 12),
##          legend.position = c(0.8, 0.25),
          strip.text = element_text(size = 8)) +
    facet_wrap(~ ID, scales = "free_y")
})
plot_PPCPK

[[1]]

. Warning: Removed 48 rows containing missing values (geom_path).
. Warning: Removed 367 rows containing missing values (geom_point).

[[2]]

. Warning: Removed 227 rows containing missing values (geom_point).

sessionInfo()
. R version 3.5.1 (2018-07-02)
. Platform: x86_64-pc-linux-gnu (64-bit)
. Running under: Ubuntu 14.04.5 LTS
. 
. Matrix products: default
. BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
. LAPACK: /usr/lib/lapack/liblapack.so.3.0
. 
. locale:
.  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
.  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
.  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
.  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
.  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
. [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
. 
. attached base packages:
. [1] parallel  grid      stats     graphics  grDevices utils     datasets 
. [8] methods   base     
. 
. other attached packages:
.  [1] mrgsolve_0.9.2     PKPDmisc_2.1.1     qapply_1.39       
.  [4] rlecuyer_0.3-4     kableExtra_1.1.0   knitr_1.25        
.  [7] gridExtra_2.3      forcats_0.4.0      stringr_1.4.0     
. [10] dplyr_0.8.3        purrr_0.3.2        readr_1.3.1       
. [13] tidyr_1.0.0        tibble_2.1.3       tidyverse_1.2.1   
. [16] bayesplot_1.7.0    rstan_2.19.2       ggplot2_3.2.1     
. [19] StanHeaders_2.19.0 metrumrg_5.57      MASS_7.3-51.1     
. [22] XML_3.98-1.20      lattice_0.20-38    reshape_0.8.8     
. 
. loaded via a namespace (and not attached):
.  [1] httr_1.4.1                jsonlite_1.6             
.  [3] viridisLite_0.3.0         modelr_0.1.5             
.  [5] assertthat_0.2.1          highr_0.8                
.  [7] stats4_3.5.1              cellranger_1.1.0         
.  [9] yaml_2.2.0                pillar_1.4.2             
. [11] backports_1.1.5           glue_1.3.1               
. [13] digest_0.6.21             RColorBrewer_1.1-2       
. [15] rvest_0.3.4               colorspace_1.4-1         
. [17] htmltools_0.4.0           plyr_1.8.4               
. [19] pkgconfig_2.0.3           broom_0.5.2              
. [21] haven_2.1.1               scales_1.0.0             
. [23] webshot_0.5.1             processx_3.4.1           
. [25] generics_0.0.2            withr_2.1.2              
. [27] lazyeval_0.2.2            cli_1.1.0                
. [29] magrittr_1.5              crayon_1.3.4             
. [31] readxl_1.3.1              evaluate_0.14            
. [33] ps_1.3.0                  nlme_3.1-137             
. [35] RcppArmadillo_0.9.800.1.0 xml2_1.2.2               
. [37] pkgbuild_1.0.6            data.table_1.12.4        
. [39] tools_3.5.1               loo_2.1.0                
. [41] prettyunits_1.0.2         hms_0.5.1                
. [43] lifecycle_0.1.0           matrixStats_0.55.0       
. [45] munsell_0.5.0             callr_3.3.2              
. [47] compiler_3.5.1            rlang_0.4.0              
. [49] ggridges_0.5.1            rstudioapi_0.10          
. [51] labeling_0.3              rmarkdown_1.16           
. [53] gtable_0.3.0              inline_0.3.15            
. [55] reshape2_1.4.3            R6_2.4.0                 
. [57] lubridate_1.7.4           zeallot_0.1.0            
. [59] stringi_1.4.3             Rcpp_1.0.2               
. [61] vctrs_0.2.0               tidyselect_0.2.5         
. [63] xfun_0.10